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A nano-shuttle consisting of two movable islands connected in series and integrated between two contacts is 
studied. We evaluate the electron transport through the system in the presence of a source-drain voltage with 
and without an rf excitation. We evaluate the response of the system in terms of the net direct current enhanced 
by the mechanical motion of the oscillators. An introduction to the charge stability diagram is given in terms 
of electrochemical potentials and mechanical displacements. The low capacitance of the islands allows the 
observation of Coulomb blockade even at room temperature. Using radio frequency excitations, the nonlinear 
dynamics of the system is studied. The oscillators can be tuned to unstable regions where mechanically assisted 
transfer of electrons can further increase the amplitude of motion, resulting of a net energy being pumped 
into the system. The resulting amplified response can be exploited to design a mechanical motion detector of 
nanoscale objects. 
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I. INTRODUCTION 

Recent experiments on coupled shuttles show intriguing ef- 
fects arising from the coupling of the electrostatics and me- 
chanical degrees of freedom (TJE). The imprints of single- 
electron effects have been observed in these systems, such as 
Coulomb blockade EH) and gate-voltage-dependent oscilla- 
tions of the conductance [5 1. The increasing relevance of nano 
electromechanical system (NEMS) is a result of their poten- 
tial industrial applications J6ll2l- NEMS offer the possibil- 
ity to realize, for instance, nanomechanical switches O, cir- 
cuits IrQ UTTI . electronic transducers lfl2l4T4ll . solar cells 031 or 
high-sensitive charge [1 1, spin 1 16 1 and mass sensors lfT71[T8l . 
as well as the general study of nonlinear dynamics of oscilla- 
tors and resonators [ 19-22 1. In particular, there is a growing 
interest in parametrically driven nonlinear systems E3l l24l . 
where a smooth change in the value of a parameter results in 
a sudden change of the response of a system. 

Our aim is to present here a theoretical study on the elec- 
tromechanics of a coupled shuttle with an external electrical 
excitation. The "smoothly" changing parameters are the in- 
tensity and the frequency of the excitation, and the response 
is the observable direct current through the system. In this 
context, we consider mechanically assisted electron transport 
through a system of movable low-capacitance nanoislands 
connected in series between two electrodes. The mechanical 
motion of the islands changes the mutual capacitance of the 
system and the tunneling processes, hence affecting the cur- 
rent through the system. We can thus say that the mechanical 
and electronic degrees of freedom are coupled. An rf exci- 
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tation allows us to study the effects of nonlinear mechanics, 
where multiple stability can be achieved by tuning the fre- 
quency and intensity of the excitation in the coupled mode 
regime. In the unstable regions the response of the system is 
greatly amplified, suggesting a practical scheme for detection 
of instabilities in the mechanical motion of nanoscale objects. 




FIG. 1: (a) SEM image of a coupled shuttle consisting of two Si-based 
nanopillars, (b) double harmonic oscillator, with coupling spring constants 
k and kg, and mechanical displacements Xl and xr, and(c) forces that gener- 
ate the displacement x in a differential element ds of the nanopillar. (d) The 
flexural mode where the center of mass is at rest. 

An example of the device that we analyze is the one devel- 
oped by Kim et al. 0, consisting of a double pillar structure 
with a gold nanoisland on top. [see the SEM image of Fig. 
[TJa)]. The device can be described in terms of two sets of 
characteristic quantities. The first one is the relative displace- 
ment of the islands, x = Xl — xr = rocoscof [see Fig. [TJb)], 
with ro being the amplitude of the oscillations and CO the vi- 
bration frequency. The mutual capacitances and resistances of 
the device depend on this quantity, affecting as well the other 
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important set of parameters: the electrostatic free energy for a 
given charge configuration, F(nij), where m, labels the charge 
state of the device. 

This work is organized as follows: in Sec. [Il] we describe 
the purely mechanical aspects of the nanopillars in a flexural 



mode. Sec. Ill then describes the electrostatics of the system, 
introducing the free energy and chemical potentials as a func- 



tion of the mechanical displacements. In Sec. IV we consider 
the coupling of the mechanical and electronic degrees of free- 
dom. A dynamic equation is derived, and we conclude with 
a master equation to describe electron transfer processes be- 
tween the contacts. We evaluate in detail the small oscillations 
limit and the shuttling regime within the Coulomb blockade 
limit. We also give an expression for the dissipated and ab- 
sorbed power by the device. Finally, we devote Sec. |V|to 
conclusions. 



II. MECHANICAL ASPECTS OF THE COUPLED 
OSCILLATOR 

The system of our interest is represented in Fig. [T] Two 
nanopillars of height h are operated in a flexural mode where 
the center of mass remains at rest in most cases, mechanically 
assisting electronic transport across the system. We consider 
first a single pillar, as in Fig. [T] (c). The beam can vibrate 
along x, with displacements x(z,t). A differential element of 
length dz and cross-sectional area A is subject to forces F x (z+ 
dz) and —F x (z) on each face, directed along x, and torques 
M y (z + dz) and —M y (z), directed along y. Balancing linear 
forces and imposing that there is no net torque l25l . we have 



F x {z + dz)-F x {z) = pAdz 



fx 
dt 2 



F x (z + dz)dz + M y (z + dz) - M y (z) = 0, 



(1) 



where M y = EI x (d 2 x/dz 2 ), E being the Young modulus and 
/, = 7tr 4 /4, the second moment of area of a cylinder. Ex- 
panding about the point z and keeping only first-order terms 
in dz, we find the Euler-Bernoulli equation (note that we are 
neglecting the damping force for now): 



EL 



fx 
dz 4 



-pA 



fx 

dt 2 ' 



with solutions x„(z,t) = [a„(cosp„z — cosh $„z)+b„ (sin P„z — 
sinhp„z)]cosco„f. With the boundary conditions x(0,/) = 
d z x(0,t) = and d 2 x(h,t) — d^x(h,t) — 0, we find numerically 
$„h = 1.875, 4.694, 7.855, . . ., with co„ = ^(EI x /pA)ffi and 
a„/b„ = -1.362,= -0.982,= -1.008,= -1.000,... l25l . 
For the first mode, we can estimate that ©o ~ 240 MHz for a 
typical nanopillar with a radius r ~ 30 nm and using a Young 
modulus of E = 150 GPa [26 1. At the metallic islands situated 
on top of the pillars (z = h), the mechanical movement thus 
can be described as a harmonic oscillator, mxi + COqX, = 0. 

We now consider the double oscillator in its coupled mode, 
as depicted in Fig. [TJb). Two harmonic oscillators of mass m 
and spring constants Icq = ma>Q, with coupling spring constant 



k = c,m(i)Q, where q is a parameter that quantifies the coupling 
of the oscillators. The dynamics of the system can be derived 
from its Lagrangian, 



L : 



1 



-mx R - -mco [x L +x R + q(x R -x L ) \ . (2) 



The problem suggests using the new coordinates, where X = 
(xl +Xg)/2 is the center of mass and x = x^—Xg, the relative 
displacement, giving 

x(t) = ro cos (cor + (p) , X(t) =Xocos(k>o? + <?'), 

where ro, Xq, <p and <p' are determined by the initial condi- 
tions. We see that the center of mass moves with the nat- 
ural frequency COo, whereas the relative coordinate moves 
with a higher frequency, CO = (1 + 2q) *J ko jm = Wkfjm, k' = 
ko(l +2q) 2 . Two limiting cases can be considered: strong 
coupling (SC), where 1, and weak coupling (WC), with 
5<1. In the SC regime, the movement of one oscillator is 
quickly transferred to the second, whereas in the WC regime, 
the movement of the first oscillator is slowly transferred to the 
second one, which is the situation we expect to encounter in 
our system. 



III. ELECTRODYNAMICS OF TWO METALLIC 
MOVABLE GRAINS BETWEEN TWO CONTACTS 

A. Free energy of movable coupled metallic grains 
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FIG. 2: A schematic picture of the double-island structure with the voltage 
sources and various capacities in the system. The two circles denote the is- 
lands with n^ R excess electrons. The distance between the islands is given 
in terms of their relative displacement, x(t). A bias V(t) is applied to the left 
contact while maintaining the right one grounded. 

The circuit diagram of our system is depicted in Fig. [2] 
two oscillating, capacitively coupled metallic islands with rii 
{hr) excess electrons in the left (right) island, resulting in a 
sequence of three tunnel junctions, i = 1, 2, 3, each charac- 
terized by a resistance and a capacitance, Rj and Q. hl,r is 
determined by the charge accumulation in the junctions, m ( - 



n L = mi -M2, 
rig = mj — rriT,. 



(3) 



A bias V(t) is applied to the left contact, while maintaining 
the right one grounded. 
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The capacitances of the (disk shaped) left and right island 
Cl,r and their mutual capacitance C2 can be expressed in 
terms of their radii, ri r, 



Cl,r — 8er L . ff ; C2 ^ 8e 



rtfR 



= C 



1 



(d—x) 1 —x/d ' 



(4) 



where d is the equilibrium distance between the islands. Here, 
e is the dielectric constant of the material surrounding the 
electrodes, e = £o£,- In our case, e r ~ 1 (air), so the dielectric 
constant is close to the vacuum one, £0. 

A full derivation of the free energy is given in Appendix [A] 
The free energy in the linear transport regime [i.e., for V(t) ~ 
0] in terms of the m, and n a given by Eq. Q reads 

1 2 1 , 

F({n a ,mi}) = -n L E C L+ ^n R E C R+n L nREcc-W({mi}); 

W({m,}) = i, { mi [V G {C Gl E CL + dlE cc )} + 
\e\ 

+ m 2 [V G (C G i (E C l - E C c) - C G 2 {Ecr - E C c)} 
+ m 3 [-V G (C G iEcc + C G2 E C R)}} , (5) 

where E C l(cr) is the charging energy of the left (right) island 
and E G c is the electrostatic coupling energy, which denotes 
the change in energy of one island when an electron is added 
to the other island. These energies can be expressed in terms 
of x, using Eq. (0), 



device oscillates between the two independent islands regime 
and a single large island regime, giving rise to a rich structure 
in the response. 

B. Transport in the Coulomb blockade limit 

We consider the electric transport within the classical 
regime studied by Kulik and Shekhter E71 . Most of the 
single-electron effects can be explained in terms of lowest- 
order perturbation theory, since higher-order tunneling pro- 
cesses, as cotunneling, are exponentially suppressed due to the 
mechanical motion of the islands. The charge state is given in 
terms of the probabilities of having ni excess electrons in the 
left island and n R excess electrons in the right island, P nL ,n R - 
The tunneling processes are described in terms of transition 
rates within the "Orthodox" model [28 1, resulting in an equa- 
tion of motion that describes the evolution of the charge with 
time. In such a picture, <?l L , m denotes the tunneling rate 
across junction j in the forward or backward direction, having 
ni and riR excess electrons in either island: 



{\-^l kBT )e 2 Ri 



(7) 



Here, R, is the resistance in the i-th junction, which depends 
exponentially on the displacement of the islands. On the fiex- 
ural mode, we have 



1 



-CL/CR 



C. 



L/R 



1 _ r JM 



Ecc — 



1 



L 2 



(6) 

We note that the charging energy of the individual islands 
do not change dramatically with the oscillations, since {d — 
x)>r L + r R implies (using r L ~ r R = r) E CL/CR < E CL / CR < 

4E® L j CR /3, with E® L j CR = e 1 j%zr L jR. On the other hand, the 
coupling energy in the regime of strong oscillations can take 
different limits: When the islands are far apart, x < and 
(d — x) 2 ^> r^rR, then E G c —> and the free energy given in 
Q is formally equivalent to the sum of the energies of two 
independent islands, 



{VgCci - n L eY 
2e 2 



Ecl- 



(VgCg2 - n R e) 
2e 2 



-E CR +f{{ mi },Vl), 



with / being a function that does not depend on x. As they 
approach, Ecc becomes larger, with Ecc ^ e 2 /12er. For (d — 
x) ~ 2r and C G \ ~ Cg2 = Cg, we separate the x-dependent 
and x-independent (g) terms, 



{4V G C G -(n L + n R )e) : 
2e 2 



Ecc + g{{mi},V%), 



which corresponds to the energy of a single island with charge 
kl + kr- Hence, in the large amplitude of oscillations limit, the 



(8) 



where is the static resistance on junction i and X is the 
phenomenological tunneling length. The electrochemical po- 
tentials or addition energies, /jp, denote the energy an elec- 
tron needs to overcome in order to tunnel across the junction 
i while keeping fixed the number of electrons in junction j 
(j ^ i), /jf (n L ,n R ) = F(mj ± 1) - F(w,-). Figure [3] shows 
schematically the six different processes across the three junc- 
tions. Using the expression for the free energy |5]), the elec- 




FIG. 3: Schematic representation of the tunneling processes in a system of 
two islands attached to the S and D leads. Six different processes are given in 
terms of chemical poten tials, f£~, i = 1, 2, 3, as indicated in the picture. 

trochemical potentials for the six possible processes of Fig. [3] 
are 
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£(n L ,n R ) 
fjf^(ni,n R ) 



l±n L )± 5z»« t ~ (C G i + C G2 S L ) ± ^ ((C t - Q) - 8 L C§) 
2 y |e| |e| 



Vff) Vr 
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CL-^CC) 



1 \ v G v(f) 

-=F«r T§««tT -tt ( c G2+C G i5«)± -j-j-CiS^ 

2 / e e 



r 



Ecr, 



(9) 



where we have defined 8 L / R = r L j R f (d — x) and made the ap- tion approach extended to multiple junctions reads 
proximation ~ r R in the second equation. The master equa- 



— ^ n L -l,n R Pn L -l,n R + ^ n L +\,n R Pn L +\,n R + ^ n L +l,n R -lPn L +\,n R -l n L -l,n R +\Pn L -l,n R + \ n L ,n R + lPn L ,n R + \ + 

+ ^nijBji— 1^V,»W— 1 — Y ^ i L ji R P"LM R - (10) 



In order to understand qualitatively the charge transport in this 
mechanically movable device, let us focus on the charge trans- 
fer from the source to the left island in the absence of a gate 
voltage and in the neutral charge state, ni = n R = 0. We can 
express the corresponding chemical potential in terms of x by 
using Q in 



Mi 



e 2 (l~2x/d) 
2C L (1+S L S R ) 



1 



Vs 

KKi-x/d) 



where we have defined the static voltage threshold \VS^\ 



\e\/2\(C L -Ci -C 2 U 5 L )| ~ H/2[C 2 U (1 - 8 L )} (note that the 
junction capacitances in the contacts are C\ = Ci - C 2 - C G i and 
C3 = Cr - C 2 - Cq2). At zero temperature, a tunneling event re- 
quires a negative chemical potential pf < 0, involving, in the 
absence of mechanical movement, a negative source voltage 
Vs < 0. For a movable system, the inequality for a tunneling 
event to occur reads 

\V s \>\Vg\(l-x/d). 



If the voltage is in magnitude just below the threshold, say 
\Vs\ = \Vfa\(l - A e ) with < A £ < 1, then a negative chemi- 
cal potential requires (1 — A £ ) > (1 — x/d), or x/d > A £ > 0, 
involving the island separating from the left contact [recall 
that x is positive when both islands separate from the con- 
tacts, as in step (2) of Fig. |TJd)] . However, the resistance in- 
creases exponentially with the distance, hence suppressing the 
tunneling process. Small oscillations are possible due to the 
elastic force, but the island will not be pushed back and forth 
by Coulomb forces. On the contrary, beyond the threshold, 
say \Vs\ = VS}(I + A E ), a negative potential no longer requires 
negative x. Tunneling of one excess electron as the island ap- 
proaches the contact then becomes possible. The direction of 
motion of the charged cluster right after the tunneling, due 



to the Coulomb forces, will be away from the contact which 
has supplied the extra electron, and thus, we may say that the 
current becomes mechanically assisted. Hence, a sharp transi- 
tion in the conductance of the system as the voltage increases 
beyond V^ 1 is to be expected. 

The same argument can be applied to the reverse jump, 
from the left island to the 'S' electrode, 



Mi 



e 2 (l -2x/d) 
2C L (l+8 L 8 R ) 



1 



V S 



K\(i-x/d) 



A negative chemical potential for the static island < re- 
quires now Vs > and |Vs| > \VS}\. If the voltage is just be- 
low the threshold, V^ 1 (1 — A £ ), a negative chemical potential 
involves as well x/d > A £ , a process then exponentially sup- 
pressed. 

It is easy to see that a similar situation occurs with the trans- 
port involving the right island and the right contact, since we 
have, by examining Eq. d9]l: 



e 2 {l-2x/d) 
2C a ,(l+5 L Sfl) 



1±- 



V s 



V^l-x/d) 



with i = 1, 3, oci = L, a 3 = R and V t g 3 = \e\8 R /2Ci. Again, 
at zero temperature, a charge transfer from the right island to 
the right contact would require the island to separate from the 
contact when the voltage is just below the threshold V^ 3 . The 

x dependency for p.^ is not as important as in the other four 
processes, but naturally, the tunneling is favored as the islands 
approach each other, x > 0. 

We compute the current in the stationary limit, i.e., when 
the transient solutions become negligible in the oscillations, 



7 DC ~ ~e Y, P «a.np 
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where the non-equilibrium probability distribution P na ,n^ is a 
stationary solution of the kinetic equation, (JTOji. To perform 
our calculations, we used materials parameters mostly based 
on typical experimental values; see Table [I] Figure [4] shows 



TABLE I: Materials parameters used in this work. (PC = personal 
communication with Prof. Dr. Robert H. Blick, AE = authors esti- 
mation ). 
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Value 


Units 


Ref. # 


Br 
HI 


6.25 




(25l 
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150 


GPa 


1261 & AF 
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PC 
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GO 


ni 
i ^ i 


"l,3 
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GO 


ni 


cS 
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aF 





C° 
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aF 
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250 


n m 


IJ 


rn 


32 


n m 


ID 


ri 


28 


n m 


ID 


cl 


45 


n m 


ID 


ro 


7.5 


n m 


AE 


L 


150 


nm 


ID 


X 


5 


nm 


AE 


m 


2xl0~ 18 


Kg 


CO CD 




250 


MHz 


CD 


£ 


10~ 2 




AE 



the numerical results of the absolute value of the normalized 



current \I^ C | using Eqs. ( 10 1 and ( 1 1 1, in units of 7 t h, the current 
at the threshold bias, Vth- The Coulomb blockade diamonds 
are apparent, a result of the discrete nature of the electronic 
charge 
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FIG. 4: Contour plot of the absolute value of the direct current |7dc as a 
function of Vs and Va- The radii of the nanoislands were set to r^ = 28 nm 
and rg = 3 1 nm. Coulomb blockade diamonds are apparent, even at room 
temperature. 



We now consider a rf signal superimposed to the dc one, 
Vs = V(t). The electron transfer between the islands to the 
right direction will occur in phase with the signal, and the 
electron transfer to the right direction between a contact and 
its nearest island will occur out of phase with the signal (see 
Fig. |5j. The sequence of electronic transport is schematically 
depicted in Fig. [5] The sign of the current is determined now 
by the initial conditions, occurring left to right in Fig. [5ja) 
and right to left in Fig. |5jb). 
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FIG. 5: Electronic transport through a double pendula structure. The direct 
current is to the right (a) [left, (b)] direction if the metallic grains approach 
the contacts (each other) when the applied bias is negative. Electron transfer 
between the islands to the right (left) direction occurs in (out of) phase with 
the signal, as depicted in bottom of (a) [top of(b)]. 

If the frequency of the rf signal is close to the natural fre- 
quency of resonance coo, the oscillatory motion can assist the 
conductance by significant displacements of the island. The 
oscillatory regime in the frame of nonlinear dynamics is stud- 
ied in the subsequent sections. 



IV. ELECTRO-MECHANICAL DYNAMICS: COUPLED 
MODE AND PARAMETRIC ELECTRONICS 

As we noted before, the capacitances and the resistances 
depend on the displacements of the islands, hence affecting 
the chemical potentials and the tunneling processes, respec- 
tively. We want to take into account explicitly the oscillatory 
movement of the islands into the electrostatics. To do so, we 
first express the free energy in terms of £q defined in Eq. (|6}: 

F(x; {mi}) = Ec L (x)T\L{{mi}) + E CR {x)T\ R ({mi}) 

+ E C c{x)r\c{{mi}), (12) 

where we have defined the set of variables r| as: 



11 L 



T\R = 



1lC = 



«L _ 
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TlR 
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n L n R 
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mj, 



[V G C G i+Vs(Cl-Ci)]- 
(VgCgi) 

'I 



7772 



VsCi] 



(Vsd)- 



Next, we expand Eci in terms of the relative displacement 
x £ i = x/d (note thatx c / < 1): 



-CL/R 



Ecc 



1 



25,. 35r + S* , 
l-8/ d+ (l-8 r )2^' 



8erf(l-8 r ; 



. 1 + 5,- l+38 r 2 
1 -5,-^ /+ (l-5 r ) 2 ^' 



(13) 
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where 8 r = rpr R /d 2 <C 1. We include the electrostatic 
force in the coupled harmonic oscillators by adding the 
term F(np,n R ;x) to the Lagrangian, L = mx 2 /2 — kx 2 /2 — 
F(ni,ri2;x). Using the leading terms in 8, for the derivative 
of Eq. (\3\ , we have: 



dF_ 
dx 




V\« y> n{n + l)x" d 1 
r R r^ (l-8,)» 



11c 



(14) 



The Lagrange's equation of motion for the relative coordinate 
x now reads: 

dF 

mx + mjx + kx = —— > mx + }nyx+ (k+A^jx = —Ap (15) 

ox 

where we have introduced the damping force, Fy — myx and 
defined 



A* 
A F 



4e</ 2 (l-8,.) 2 
e 2 

8erf(l-8 r ) 



35,. 



28 r 



rt 

fL 



Equation ( p~5] > is known as a damped Mathieu equation PP . 
It includes a standard harmonic oscillator driving term, and a 
parametric modulation term, which is a variable spring con- 
stant, Ajfc. We note that A* and Af depend in general on time, 
as the shuttle oscillates causing a charge transfer. A rough es- 
timation for a realistic system gives Ap ~ 0.2 — 1 pN. T\l,rc 
depend on the number of excess electrons on each of the is- 
lands, n L f R . As we will see below, in the shuttling regime, 
n L i R is a function that oscillates with time, as the mechanical 
movement of the islands assist the electronic transport through 
the device. Thus, A^ and Ap are time-dependent functions, 



and Eq. ( 15 1 is a modified Mathieu equation, which can be 
treated numerically |[32l . In the following subsections, we 
consider two different limits to understand the dynamics of 
the system. 



A. Small oscillations limit in the linear regime 

We consider first the small oscillations within the classical 
circuit limit. In the linear adiabatic regime, the charge bal- 
ance in the islands follows the excitation given by the applied 
voltage, V(t). From classical circuit theory, we have that the 
applied voltage equals to the sum of the voltages that drop on 
each junction, and the net current through each of the junc- 
tions is the same, 



Vfc> = V(t) = £ 



9£. 

Ci 



RiQ 



1} 



(16) 



We can solve the above system of equations for = — |e|m,- 
and get the charge on each island, Qi — — \e\nf. Qp = q\ — q2\ 
Qr = 12 = <?3- When the flexural modes of the nanopillars are 



excited, the mutual capacitance and the resistances become 
sensitive to the displacements, xp and x R . In a flexural mode 
in which the center of mass is at rest (Xq = 0), as in Fig. |TJd), 
the resistances are given by (fij). The mutual capacitance C2 
determining the coupling of the metallic islands depends as 
well in their separation, whereas C13 can be considered as 
constant, 



C1.3 =2 C 



1,3' 



C 2 (x) 



L 2 



1 —x/d 



Further, we will make the approximation that R®Cf is a con- 
stant, but the resistance R2 is twice the resistances Ri 3. We 



set R 



■Ri 



: R\/2 and C\ = q 



2C\ 



C, consistent with 



previous results |3|. From now on, we express the relative 



coordinate in units of X, x = x/X. Using ( 16 1 we get 



02 



CV(t) 



41 



2(\-x/d)(\+e ix l 2 Y 
giving a compact expression for Qp 



CV(t)e 3x / 2 
2(1 +«a*/2) : 



-Qr, 



Ql 



CV 



u 3x 

tanh 



x 3x/4 
d 



3CV 



2 2 3 3 

-XT X 3 

d 16 



(17) 

At lowest order in x, we have a symmetric system whose is- 
lands have no net charge at x = 0. Then, for x < (so the pil- 
lars are away from each other, getting close to the contacts), 
a net charge with negative sign is induced in island L and a 
positive one in R. Note that the potential is, by convention, 
negative on the left and positive on the right contact. Revers- 
ing the potential V will cause a change of sign in Qp and Q R , 
as expected. The two nonlinear terms on the right break the 
left-right symmetry of the system. 

The dynamics of the nanopillars will consist of a set of os- 
cillators experiencing the electric field, V(t)/L. We follow the 
approach given by Ahn et al. Il22l to investigate the electrody- 
namics of the system. We find the equations of motion of the 
relative coordinate x by setting m = m R ~ mp and substituting 



Eq. (17 1 into Eq. (15 1 



Y • -2 
— x+x — —a sin" cot 

COo 



2X 2 3 3 
X -Y X "16* 



Here a = 3CVq/4L;«A-C0q, a dimensionless forcing parameter 
which account for the ratio of the electric (F e ~ CVq/L) and 
mechanical forces (F m ~ mXcOg = kX). We note that for a typ- 
ical Si nanopillar, F m - 50 - 60 pN and F e ~ 1 - 5pN. We 
have also rescaled the time, X = (Do? . CO is the frequency of the 
excitation, expressed in units of COo. This is a non-linear equa- 
tion that corresponds to a forced and damped oscillator, where 
the forcing terms depend on the coordinate itself. At first or- 
der in x, we obtain a modified Mathieu equation, which gives 
instability regions when the excitation is strong enough. We 
are, however, interested also in the weak excitation regime, 
in which the non-linear terms and non-linear effects such as 
Coulomb blockade could play a critical role. 

Following the Poincare-Lindstedt method, we parametrize 
the damping and the forcing using a small arbitrary e (£ <C 1), 
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Y~eyi, oc~eoci, 



x + x + e yii + oci sin" cot 



2X 2 3 1 

x x x 

d 16 



= 0. (18) 



"slow" time T|. At first order in e we get: 



d 2 Xl 

IF 



d 2 XQ 



dx 



Xl = -2^-25^-Yt^ 



Thus, we can consider two different time scales, the 
"stretched" time, z = COT, and the "slow" time, r| = £T. The 
time derivatives are now expressed in terms of these new times 
as 



cxi(l — cos2z) 



xo - 



2Xxq 



H 

16 



(23) 



> d 2 x _ d 2 x 



d 2 x 



dx dx 

i = co^-+e^— ; x = oo ^+2eco^-^- + e"^r. 
6z or) dz z or|dz dr| z 

We expand x in terms of e, 

x(z,T|) ~Xo + £X\+... 



(19) 



(20) 



Without loss of generality, we may ask that x\ satisfies x\ + 
x\ =0. We substitute the general solution ( p2] i for p = 1 into 
the above expression and arrange terms in sinz and cosz (see 
Appendix B 1 1 to get: 

.dA 



B- 3 ^B(3A 2 + 5B 2 ); 
64 v ' 



and, likewise, seek for solutions that correspond to harmonics 
of the natural frequency C0q: 



2— = -fiB- 
dr) 



26a,- 



3aj 

"64 



A(A 2 + 3B 2 ). 



(24) 



co ~ p + eS a 



(21) Note that equilibrium points of ( 24 1 correspond to periodic so 



lutions of our forced oscillator, and the norm of the solutions 



where p is an integer or fractional number and, finally, substi- is conserved, i.e.: A +B = x\ 



tute (19 1, (20 1, and (21 



into 



18 1, neglecting terms of (9(e ), 



which gives, after collecting terms at lowest order in e, 

1 d 2 XQ 



If a dc signal is superimposed, V(t) — Vb(sinC0T + p) with 
P = V c ic/Vo, Eq. ( 18 1 will now read 



p 2 dz 2 

giving a general solution for xq, 



x =0 



8 <j Yii + ajsin C0T + 2psincoT + p ] 
jc = 0. 



2X 2 3 -i 
X -V X ~16 X 



x (z,r|) = A(r|)cos - +B(r))sin-^ 



(22) 



The constants of integration, A and B are functions of the 



Proceeding in the same manner as before, we find six extra 
terms (terms in P or p 2 ). We consider the stretched and slow 
time, and expand x and co in terms of e, to get, for the p = 1 
case (see Appendix |B 2\: 



dA 
dr| 
dfl 
dr| 



-YiA 
-Yi«+(6 



a. 



^ {A 2 + 3B 2 } 2 J^ B{A 2 +B 2 ) ^ B{M 2 + 5B 2 } 
a 64 lzo 



«ip 2 ) A 



2aipA, 



AB 



27a ^ 2 A(A 2 +B 2 ) + 3a ±A(A 2 + 3B 2 ) 
64 V ' 128 V ' 



(25) 



and for the p = "any" case, (see Appendix 



B3i: 
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,dA 
dr) 



,dfl 
dr) 



-A(y 1 -pa 1 5,, 2 )- J B^-^(2 + 5 p , 1 +4p 2 ) 

3ai 
~256 



AB 



— S p , 3/2 + (A 



[fi(A 2 + fi 2 )(6 + 12p 2 + 38 p ,i) + fi(3A 2 -fi 2 )8 p , 2 +6A(A 2 +fi 2 )8 p , 2 +2A(A 2 - 3B 2 )S P>4 ] 



-B(y 1 + Pa 1 8,, 2 )+A^-^(2-8 p , 1 +4p 2 ))+(A 2 -B 2 )^8 p3 / 2 -AB^8 p ,3 
~ [A(A 2 +B 2 )(6+ 12p 2 - 38 p .i) - A (A 2 - 3B 2 )8 P , 2 +6B(A 2 +B 2 )8 P , 2 +2B(3A 2 -B 2 )8 P , 4 ] 



(26) 



The stability of the solutions of the equation above can be 
investigated using numerical methods l32l . 

In order to understand qualitatively the electromechanical 
motion scenario, we consider, first, the case p = 1 and P = 0, 
where Eq. (25 i reads (A = dA/drj) 



2A = -YiA-(y 2 -Y3)B- 
IB = -YiB + (Y2+y 3 )A- 



-y 4 fl(3A 2 
-Y4A(A 2 4 



h5B 2 ) 
3fi 2 ) 



where we have defined y 2 = 28(o — oci, y? = ai /2, and y 4 = 
3ai/64. It is easy to see that the transition curves for the sta- 
bility of the trivial solution are y 2 = ±y;, or ai = 48a,, 48 m /3. 
Along this curves (broken lines of Fig. |6| the stability of ro 
= changes. To gain in simplicity, we transform to polar 
coordinates in the A-B phase plane, by setting A = ro cos cp 
and B = rosincp. The amplitude r 2 , = A 2 +B 2 and the phase 
cp = arctanfi/A now satisfies: 



2r = -yir + r sin2(p(y 3 -y 4 rg) 
2<p = y 2 +y3cos2(p+y 4 ro(3 — 2cos2(p) 



(27) 



We seek equilibria of the "slow flow" ( |27"| ). A solution in 
which ro and cp are constant represents a periodic motion of 
the nonlinear Mathieu equation, which has the frequency of 
the forcing function. Such equilibria satisfy ro = cp = 0. Ignor- 
ing the trivial solution ro = 0, the first equation of ( |27] > with ro 
= requires 



sin2(p ; 



Yi 



Y3-Y4^o 



In the absence of damping yi ~ 0, find equilibria at cp = 0, 7c/2, 
71, 3tc/2. The second equation of (27 1 with cp = then implies 



y 2 +Y3Cos2(p 



Y2TY3 



y 4 (3-2cos2cp y 4 (3=p2)' 



For a nontrivial real solution, > 0. In the case of cp = 
or Jt, cos 2cp = 1 and nontrivial equilibria require — ys — y 2 > 
or 48(o < oti . On the other hand, for cp = jc/2 or 37C/2, 
cos2cp = — 1 and nontrivial equilibria require Y3 — y 2 > or 
48(o < 3oti. Since b a = U\/4 and 5 a = 3oti/4 correspond to 
transition curves for the stability of the trivial solution, bifur- 
cations occur as we cross the transition curves in the 8(o-OCi 
plane (see Fig. [6]). Keeping ai fixed as we shift 8 a from the 



right across the right transition curve, the trivial solution ro = 
becomes unstable and simultaneously two branches of sta- 
ble solutions are born, one with (p ~ tc/2 and the other with 
9 ~ 3tc/2. This motion grows in amplitude as 8(o continues to 
decrease. When the left transition curve is crossed, the trivial 
solution becomes stable again. This scenario can be pictured 
as involving two pitchfork bifurcations, as depicted in Fig. [6] 



a, 




"<p ~ 37C/2 



FIG. 6: Schematic representation of the p = 1 tongue in the parameter space 
spanned by S ra and a\. Between the transition curves (broken lines), the 
trivial solution is unstable and two stable solutions are born, (p = 0, % (solid 
curves). The insets depict phase portraits in the A-B plane: below the lower 
transition curve, only the trivial solution exists. In between both transition 
curves, the trivial solution is unstable, and another two solutions are bom. 
Above the second transition curve, ro = is stable again, and the other two 
solutions become unstable. 

A finite damping yi shifts slightly the position of the stable 
equilibria in the A-B plane, which become stable spirals, as 
represented in the insets of Fig. [6] The damping also "shrinks" 
the region of instability (shaded area of Fig. 15), lifting it away 
from the origin in the parameter space. 

In the classical limit, we can have an idea of the resulting 
direct current through the system. The time-average direct 
current is obtained by integrating over a period the current 
across one of the three junctions of Fig. [2] 



co f'o+ T V(t)e x 
A%Rj t0 l+ e 3.v/2- 



(28) 



We study ( 28 1 in terms of the coefficients ro and cp, x — 
ro cos (cof —9). We find that the absolute value of the direct 
current reaches a maximum when cp = 0, K and ro = 2. 
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Figure[7]shows numerical results of I<x c in the p — 1 tongue, 
where 0Ci is in log scale. Inside the tongue, the trivial solution 
is unstable, and two stable solutions appear at P % /2 and P^ % /2 
in the A-B plane (see inset). The trajectories stay either in the 
upper or lower semiplane; thus, the initial conditions deter- 
mine the point of stability: If (p(0) > (< 0), then the phase 
portrait in the A-B plane reaches P n / 2 (P371/2) an d the electron 
transport is right to left (left to right), using the convention of 
Fig. [5] At these points, the nanopillars are oscillating with the 
natural frequency of the oscillations, mechanically assisting 
the electronic transport. 



the dynamic equation, [Eq. [29) in the stationary regime. In 
the low bias limit, n(t) is a sinusoidal function that follows 
the excitation V(t) (thin solid curve of Fig. The charge 
has only a small probability to be transferred across the de- 
vice, following the bias. However, after some critical bias, 
the charge transfer occurs mostly at the points of maximal de- 
flection, commonly termed as shuttling regime. We obtain 
numerically n(t), resulting a square wave correlated with the 
excitation V (t ) (thick solid curve of Fig. [sj, 

n(f) ~n av +4«o(cosc0of — cosScOoO/ 71 ; (30) 



B. Oscillations in the shuttling regime 



So far we have studied the limit when the charge in the 



metallic islands is given by Eq. ( 17 1, i.e., the charge on the 
islands Ql,r changes continuously, following the excitation 
V(t). As the size of the metallic islands shrinks, however, 
we may reach the discrete limit, where single-electron ef- 
fects such as Coulomb blockade become important. In the 
Coulomb blockade limit, as we have seen in Sec. |HIB[ an ex- 
tra electron can only be added to the island if enough energy 
is provided by the external sources to overcome the Coulomb 
repulsion between the electrons. The equation of motion for 
the relative coordinate reads: 



-yx- 



eV(t) 
klX 



n(t), 



(29) 



where n(t) =ni(t)— n«(f). 

The evolution of the probabilities of having Hl excess elec- 
trons on the left island and tig on the right island is given by 
Eq. ( 10 1 , which is solved by direct integration. Figure [8] 
shows numerical results of the function n(t) obtained solv- 
ing simultaneously the master equation [Eqs. f7j)-([T0|)] and 



where n av and no are obtained after averaging over a large 
number of simulations and depend on the input parameters. 
We note that the charge transfer occurs in the points of maxi- 
mal deflection during an effective contact time, in accordance 
with Weiss et al. [29). 

Inserting the expression ( 30 1 into p9| ), we aim, as before, 
for oscillatory solutions, xo ~ A(r|)cosz/p + B(r[)smz/p, 
bearing in mind the following linearized equations for the co- 



efficients A and B (see Appendix B 4 1 



dA 28(n , „ a', . „ . 

2— = -y 1 A-—B+n sv a l S Pt i+n -±{28 Pi 2-8p,4) 



if. = _ YlB + ^A+aip«o, 
dr| p 



(31) 



where we have defined a' = eVo/LkX and a'j = ea'. As be- 
fore, a' can be viewed as the ratio between the electrical and 
mechanical forces. We find the equilibrium or stable points by 
solving Eq. ( 3 1 1 with A = B = 0. In the absence of Vd c ((3 = 0), 
non-trivial stable points are found only around subharmonics 
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with p = 1, 2, or 4, 

04eq,5 e q)p=l 
(A eq ,-S e q)p=2 
(A eqi B eq ) p= 4 



(Yi,25 m ) 



(Yi,8a 



a'«o 

3(y 2 + §») 2 
a'no 
6(7? + (8,0/2)2) 



Yi, 



(32) 



In contrast, under a finite dc bias, nontrivial stable points are 
found for any frequency at 



(A e q,-Beq)p : 



apno 



' 28(o 

v2_l COS M2l - ' I 1 



(33) 



y 2 + (28a,/ P ) 2 V P 
We stress that this result is valid in the limit of large oscil- 



lations. To estimate the range of validity of Eqs. (32i and 
( |33| l, we consider the limit of small oscillations, for which the 
charge on the islands "follows" the mechanical motion, with 
an amplitude proportional to the amplitude of the oscillations, 

r = VA 2 +B 2 , 

n(t) ~ a„rosinCDo?. 



In this limit, Eq. ( 3 1 1, in polar coordinates and for p ^ and 
p ^ 2 reads 

r\ = -^-rl + a[$a n Ar , 

We find a change in the stability of ro = when the forcing 
parameter exceeds a threshold, a'jpa,, > a[ h , with 



ro changes from a stable to an unstable spiral as the forcing is 
increased beyond oc( h /(pa„). In this situation, the amplitude 
of the oscillations could reach the Fowler-Nordheim tunneling 
limit ll33l . with a subsequent enhancement of the direct cur- 
rent due to field emission, marking the limit of validity of the 
present model. 



finite dc bias cases, the amount of energy pumped into the 
system per cycle in steady state [at (A eq , B eq ) p ] is 



<W a > 



root' 
12 



[6pnocoscp + (6n a v8 p ,i +2no8p,2-«o8 p ,4)sin(p] . 



If this amount is larger than the dissipated power, (W a ) > 
(Wdi s ), self-sustained oscillations are expected. This may 
occur when max{oc'n av ,o, oc'«op} > Y r 0> with the appropriate 
phase (p. In particular, for the branch with <p ~ 0, the condition 
reads oc'p«o > Y r o- In other words, self-sustained oscillations 
may occur for a wide range of frequencies if the appropriate 
values of a and p are met. 



V. CONCLUSIONS 

We have theoretically studied a coupled shuttle consisting 
of two oscillating nanoislands connected in series between 
two contacts. We express the chemical potentials in terms 
of the relative distance of the islands, /jp [x(t)], and numer- 
ically integrate the master equation to obtain the direct cur- 
rent through the system. Under a dc bias, Coulomb block- 
ade diamonds were obtained. Adding an rf signal, we ana- 
lyze the response within the context of nonlinear dynamics. 
We study qualitatively and quantitatively the structure of the 
mode-locked tongues in the parameter space. Parametric in- 
stabilities are observed in a range of applied voltages and fre- 
quencies, where resulting small mechanical oscillations are 
amplified. In this instability region, an rf signal can be ex- 
ploited to parametrically amplify the response to a gate exci- 
tation. Hence, we propose a practical scheme for direct de- 
tection of instabilities in the mechanical motion of nanoscale 
objects. 



Acknowledgments 

We are grateful to R. H. Blick and C. Kim for enlightening 
discussions. This work was supported by the Spanish Ministry 
of Education, program SB2009-0071. 



C. Dissipated and absorbed power by the oscillators 



We focus now on the dissipated and absorbed power by the 



shuttles. According to Eq. (29 1, the (unitless) power loss per 



unit cycle of duration T will be given by 

(W dis )=j(x 2 )= 1 -jr 2 . 

Similarly, we can obtain the absorbed power per unit cycle by 
averaging over a period the pumped energy of the electrostatic 
force given by the last term of Eq. ([29] 



(W a ) ~ (a'n(f)x(f)(sinO)f + p)). 

Due to the correlation between charge fluctuations n(t) and 
the velocity of the nanopillars, x(t), a positive amount of en- 
ergy may be pumped into the system. For instance, for the 
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Appendix A: Free energy of a double coupled metallic island 

In this appendix we derive the free energy of a double 
metallic grain system. In general, for a system of N con- 
ductors, the total charge on each node j is the sum of the 
charges on all of the capacitors connected to node j, —enj = 
~ e Hk m k = HkCk(Vj — Vk), where Vj is the electrostatic po- 
tential of node j and ground is defined to be at zero potential. 
The charges on the nodes are linear functions of the potential 
of the nodes, Q = CV, where C is the capacitance matrix. For 













R2C2 






V, m 


V 2 m 2 




FIG. 9: A schematic picture of the double-island structure with the voltage 
sources and various capacities in the system. The three junctions are charac- 
terized by nij and Vj, j = 1, 2, 3. The islands are coupled to each other, with 
a mutual capacitance C2, as well as to a gate voltage, Va with capacitances 
Cai and Cq2, respectively. Also, they are coupled to the source-drain leads, 
with capacitances C\ and C3. 



the system depicted in Fig. [9] we have: 



Ql 
Qr 



c 2 (v L - 
c 2 (v R - 



Vr) 

■v L y 



-Ci(Vl- 

-c 3 (v r 



■Vs)- 

v D ) 



c Gl (v L - 

-Cgz(Vr - 



Vg) 

-V G ),(A1) 



where is the "residual" charge in dot i when all potentials 
are grounded. We can write this in the form Q + a = CV, with 
C being the capacitance matrix, 



QL + CiVs + CciVc \ 
Qr + C 3 V d + C g2 V g ) 



c L -c 2 \ I v L 
-c 2 C R \v r 



where Cl = C\ + C 2 + Cq\ and Cr = C3 + C 2 + Cq 2 . We can 
thus set our node voltages in terms of the capacitances, 




1 



CiCr — C\ 



C R C 2 \( Q L +CiVs+C Gl V G \ 
C 2 C L j \Q R + C 3 V D + C G2 V G ) 



The electrostatic energy of the system can now be calculated, 
U = VCV /2, with V S = V D = V G = 0, 



U(n L ,n R ) = 



ClCr 



C 2 



1 



1 



To calculate the free energy, we first calculate the work per- 
formed by the external sources to achieve a configuration with 
»i = mi — m 2 and iir = m 2 — m 3 electrons in the system. We 



consider Eq. ( Al 1 in terms of the voltages on the junctions, Vj, 
j = 1, 2, 3. For simplicity, we consider, first, the case without 



Vg, le., C G1/2 = 0: 
1 



Vi 



[C2C3V (t) + em\ (C 2 + C3) — em 2 C 3 — em 3 C 2 ] 



V 2 = - [CiC3V(t) -em\Cs -em 2 {C^ +Ci) -em-sC{\ 



V3 



[C\C 2 V {t) - em\C 2 - em 2 C\ +eni3(Ci +C 2 )} , 



(A2) 



where we have defined E = ClCr — C\. We now calculate the 
work done by the external source, V(t) = Vs. For instance, 
the work done by Vs for one electron to tunnel through the 
first junction, m\ — > m\ + 1, is given by 8W1 = V$qi, where 
8q\ = — \e\ +C18V1. From Eq. (A2i, we see that 



SVi = \e\ C ^- 8W1 = -\e\V s ^. 



In general, the work performed by the external source V$ for a 
tunneling event through the /-th junction is 



W i = -\e\V i 



{Ziik) 2 CjCk 



2E 



where e,-;j is the Levi-Civita antisymmetric tensor. For the 
finite C G \j 2 case, we have (now Cl = C\ +C 2 + C G \ and Cr = 
C 2 +C 3 +C G2 ): 

Vi = -[CR{C L -C l )V(t)-(C R C Gl +C 2 C G2 )V G + 

G 

+ \e\(miC R -m 2 (Cg-C2)-m 3 C 2 )]. 

We can now write an expression for SW, in terms of the capac- 
itances, 



?\Vs 



8W 2 = -\e\Vs 



C 2 Ct,+C G \Cr + C g2 C 2 



C\Ct,+C\C g2 



8W3 



?\Vs 



CiC 2 



Likewise, the work performed by the gate in the tunneling 
process through the /-th junction, 8Wg/ = V G [5q G \ {nit — > m,- + 
1) +8q G2 (mi -> trii + l)]: 

C G iCr+C g2 C 2 



5W G i 
5W G2 



Wg 
\e\V ( 



C G2 (Cl — C 2 )— C G \ (Cr — C 2 ) 



W G3 = -\e\V G 



C g2 Cl + C 2 C G \ 



(A3) 



The electrostatic free energy is then given by F = U — W, with 
W({nii}) — £m,(8W,' + 8W&), and, bearing in mind Eq. 
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2 ( 2 2 

F({n a ,m,}) = 6 I C R ^+C L ^-+C 2 n L n R - ^ [V G (C G iC* + C G2 C 2 ) - V S [C 2 C 3 + C G iC« +Cg 2 C 2 ]] 
C l C r -C$ { 2 2 \e\ 

- ^ [v G (c G i (c* - c 2 ) - c G2 (c L - c 2 )) + y 5 CiC 3 + CiC G2 ] + ^ [v G (CiC G2 + c 2 c G o + v s c,c 2 ] 



(A4) 



Appendix B: Transient equations 
1. rf excitation 



We want to arrive from Eq. (23 1 to Eq. (24 1. We use xq 
A(r|)cosz + 5(r|)sinz andEq. (19 : 



iXo 
dz 



= — Asinz + Bcosz; 



dz 2 



= — A cosz — fisinz; 



d 2 x dA dB 

= -^sinz+ ^cosz. 
dzdr| dr| dr) 



Also, we need the terms in xq, Xq, and Xq all multiplied by 
(1 — cos2z). For the linear terms, we have to evaluate (1 — 
cos2z)xo, for the quadratic terms, (1 — cos2z)jcq, and, finally, 
for the cubic ones, ( 1 — cos 2z)xq, 



(1 — cos 2z) (A sin z + B cos z) 
(1 — cos2z)(Asinz + Bcosz) 2 
(1 — cos2z)(Asinz + Bcosz) 3 



A B 

— (cos z — cos 3z) + -(3sinz — sin3z) 

A, . AB r . B 2 , 

— (1 — cos4z) + — (2sin2z — sin4z) + — (3 — 2cos2z + cos4z) 



— cosz cos3z cos5z 

4 V 2 2 



3AB 2 / 3 „ 1 

cosz cos3z+ -cos5z 

4 V 2 2 



3A 2 B / . 1 . „ 1 . ' 
— , — sinz+ -sin3z — -sin5z 
4 V 2 2 

B 3 /. 5 . , 1 . c 
- — 5sinz- -sin3z+ -sin5z 
4 V 2 2 



r 



Now all we need to do is to substitute the above equations into p, P 2 ), 

Eq. ( f23| l and arrange terms, i.e., we obtain an equation of the 

form: 



'dz 2 



+X\ — (...) sinz +(... ) cosz + nonresonant terms. 



For non-resonant terms, we require the coefficients of sinz and 



cosz to vanish. We then get Eq. (24i 



2. Superimposed DC excitation 



We consider Eq. (25 i and as before, expand the solutions 



in terms of e, using Eqs. (19i, (20i, and (22i and arrive at 
a similar set of equations, only now we have extra terms (in 



d 2 ;q 

dz 2 



28, 



d 2 * o d 2 x 
dz z dzdr) 



1 dxo 

e dz" 



2a(sin z z + 2psinz + p 2 ) x Q 



21 



-x - 



3 x 3 
A () 



16 



= 
(Bl) 
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The last term is a tedious one, giving nine terms, six of which 
are new. We use 



3. Subharmonics in the general case 



2a sin zxq 



a 



[Acosz + 3fisinz] 



4A,a 



2 2 

-sin zx 



d 

— sinVo ^ -— A(A 2 + 3B 2 )cosz- 
16 64 L v y 

+fi(3A 2 + 5fi 2 )sinz] 

4apsinzxo ~ 2apB 

d 

9ap 

16 u ~ 

2ap 2 x ~ 2ap 2 (Acosz + fisinz) 
4A,ap 2 2 2Xap 2 



8A,ap . , 

— sinzx 

d 

12a . , 

■ smzx 



B[(A 2 + 3B 2 ) sinz + 2ABcosz] 
B(A 2 +B 2 ) 



d 

3ap 2 

j 

16 



d 
9ap 2 
32 



-(A z +fi 2 ) 
(A 2 +5 2 )(Acosz + 5sinz) 



We consider the general case in which the driving frequency 
is given by CO ~ p + e8 m and x ~ xq + ex\ , giving: 



i(f) 
x(t) 



dx dx dx : 

l-e— +ep— + 0(e' 

^ dz 



,2\ 



+e5(o— +e— 
dz dz dr) 

2 d 2 ^ d 2 X d 2 X 2 d2jf l , rM 2\ 

P d? +2e5 ^ d? +2e ^ + ep d^ + 0(£ } 



(B3) 



So now Eq. ( 25 i is, to lowest order in e, 



Substituting into Eq. (Bl I and again arranging terms in sinz 
and cos z, 



2^ = - Y A-('28 a) -^ t -2a^)6> 



2]d 2^P_ (a2 + 3z}2) 



3a 
64 



B[3A 2 (l+2p 2 )+B 2 (5 + 6p 2 )] 



~dB / „ a „ „,\ A 4A,ap , 

2 a, = -yB + {2K- 2 -2^ 2 ) A+ ^AB- 



3a 
64 



A[A 2 (l+6p z )+3B 2 (l+2p 2 )] 



2 d2 *0 , n 



z z 
jt = A(r))cos - +B(r|)sin-. 

P P 



(B2) 



Substituting the expression for xq into Eq. ( 25 i to first order 
in e, we get 



25 m 
P 



z z 
A(r))cos - +5(r|) sin - 

P P 



dA(r|) . z dB(ri) z 

sin cos — 

dr) p dr) /? 



1 

Q 



z z 
A(r|)sin - +B(r|)cos - 

P P 



2a(sin 2 z + 2Psinz + P 2 ) 



A(r))cos - +B(r))sin ■ 



3 

16 



z z 
A(r|)cos - +Z?(r|)sin - 

P P 



1 3 X 



= 



~d 



z z 
A(r|)cos - +B(r\) sin - 



(B4) 



Next, we need to evaluate the last nine terms of the equation 
above, which involves trigonometric operations. The relevant 
contributions for the dynamics of the system are the terms pro- 
portional to sinz/ p and cosz/p (resonant terms). These are: 



. 9 a 
axo sin z = — 



z z 
A(2 + S„ i)cos — 1-5(2 — 8„ i) sin — 

P ' P 



Pa 2 xo = pa 2 



z z 
A sin — hBcos — 



A,0Ua 9 A-a 

^r slirz = i6^ 



(B 2 -A 2 )cos-+ABsin- 
P P 



2paxosinz = Pa 



A sin — hBcos — 
P P 



Jp,2 



A,ap , . A,ap 

— -x sinz = — — 
d Ad 



2Aficos- + (A 2 -B 2 )sin 
P P 



Jp.3 



14 



Xap 2 2 



2d 



Xq sinz — > no contributions 



16 - 



9ap 2 
~~64~ 



(A 



2 , B 2, 



A sin — hBcos — 
P P 



^sin Z 



3a 
128 



+fisin (^3(A 2 +B z )(l - ^.i) - -(3A 2 -fl 2 )S p , 2 



A cos i hi(A 2 + fi 2 )(l - |8 Pl i) - ^(3A 2 -fi 2 )S p , 2 



3 ax 



16 



sin 2 z 



3a 
128 



A sin i (3(A 2 + fi 2 )8 p , 2 + (A 2 - 3fi 2 )S p . 4 ) 



-ficos - (3 (A 2 + fi 2 )8 p . 2 + (3A 2 -B 2 )8 P , 4 ) 



As we did before, we can now arrange all the terms in Eq. 
(B4i as coefficients of sinz/p and cosz/p, resulting in our 



desired equation, 



2 3r| 



A (y _ p a8; , : . - S ( ^- "(2 + 8,,! +4P 2 )) +AB^8 / ,,3/2 + (A 2 -B 2 )^8 p ,3 



3a 
256 



[B(A 2 +B 2 )(6 + 12p 2 + 38 p ,i) +B(3A 2 -B 2 )8 P , 2 + 6A(A 2 + fi 2 )8 p , 2 +2A(A 2 - 3fi 2 )S p , 4 ] 



,3B 
'3r, 



= -fi(y+pa8 p , 2 )+A 



28 a 



a 



--(2-8 p . 1+ 4p 2 ) \+{A 2 -B 



hx 
16d' 



V3/2 



-AB— 8,3 



3a 
256 



[A(A + Zr) (6 + 12p 2 - 36,,!) — A(A — 3fi 2 )8 p , 2 + 6B(A 2 +5 2 )8 P , 2 + 2B(3A 2 - fi 2 )S Pj4 ] 



(B5) 



4. Transients in the Coulomb blockade limit 

We substitute again the general solution, x(t) ~ xq + Exi, in 
Eq. (29 1, along with the expression for (n(t)) given in (30i. 
The equation of motion now reads 



As before, we evaluate the equation of motion with Xo(t) ~ 
A(r|) cosz/p + B(y\) sinz/p, using (Bli. Inserting this into 



Eq. ( |B6[ ) and considering the resonant terms, with CO/ = z and 
©of = z/p, we obtain 



a 2 . 



d 2 x Q 



d 2 xo 1 3xo 



~^r~^ • 2. . — • 28(„ — t^~; 

dz 2 dzdr| dz 2 Q\ az 

a'j (sincor + P)(n av +4n (cos (Oof - cos 3coo? )/%) 





(B6) 
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o §(0 . z 
2 — A cos — 

P P 



O 8(o„ . z 
■2 — Bsin — 

P P 



dA . z 
■ 2 — sin — 
dr) p 



z , . z z 

-2 — cos yiAsin — hyiocos — 

dr| p p p 



-a[ p«o cos - 



■a\n 



l«av 



sin— S„ i 
P 



1 , . z a 
-a^osin-bp.2 

2 p 



1 , • Z/ S 
-a 1 « sin-(6p.4 
6 p 



a'j p« av 

'§/,2) 



(B7) 



AiTanging the terms in sinz/p and cosz/p, we get Eq. (31 1. 

In the small oscillations limit, n(t) "follows" linearly the 
mechanical oscillations of the islands, n{t) ~ «o sin ©of, giving 
a different set of equations, 

2^ = -YiA-^B + cciPno, 



dr) p 

dB 28ft, a', „ 

2— = - yi B+^A-«o^8 p ,2. (B8) 
dr| pi' 
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